C	Copyright Eric Fuller, 2001-2007
C	This program was derived from an example in
C	Numerical Recipes in C. Eric Fuller extended
C	the example program to 2 dimensions.

	PROGRAM PFIT

	INCLUDE 'pfit.cmn'



	REAL*8 X1(NMAX),Y1(NMAX),SIG1(NMAX)
	REAL*8 X2(NMAX),Y2(NMAX),SIG2(NMAX)
	REAL*8 A(MMAX),SIGA(NMAX),Z1(NMAX),Z2(NMAX)
	INTEGER*4 NDAT,MA,AI(MMAX)

	REAL*8 COVAR(MMAX,MMAX),BETA(MMAX),AFUNC1(MMAX),AFUNC2(MMAX)

	REAL*8 CHISQ,PHI,RHO
	INTEGER*4 I,J,K,L,M,MFIT

	EXTERNAL FUNCS
	EXTERNAL GETPARAM,GETDATA,MODEL,CALCCHI,STATS
	EXTERNAL SORT,BUILD



	MA=18

	CALL GETPARAM(MA,A,AI)
	CALL SORT(MFIT,MA,AI,COVAR,BETA)
	CALL GETDATA(PHI,RHO,X1,X2,Y1,Y2,SIG1,SIG2,NDAT)



	CALL BUILD(NDAT,MFIT,X1,X2,AFUNC1,AFUNC2,MA,PHI,
     *  Y1,Y2,SIG1,SIG2,
     *  COVAR,BETA,A,AI)



C	SOLVE MATRIX

	CALL GAUSSJ(COVAR,MFIT,MMAX,BETA)

	J=0
	DO L=1,MA
	  IF (AI(L).NE.0) THEN
	    J=J+1
	    A(L)=BETA(J)
	    ENDIF
	  ENDDO

C	RECOVER THE PARAMETERS AND CALCULATE THEIR UNCERTAINTIES

	CALL COVSRT(COVAR,MMAX,MA,AI,MFIT)

	DO I=1,MA
	  SIGA(I)=SQRT(ABS(COVAR(I,I)))
	  ENDDO



	CALL CALCCHI(CHISQ,NDAT,
     *  X1,X2,AFUNC1,AFUNC2,MA,PHI,A,
     *  Z1,Z2,Y1,Y2,SIG1,SIG2)



	CALL STATS(NDAT,MFIT,X1,X2,Y1,Y2,Z1,Z2)



C	SHOW RESULTS

	WRITE(*,53)

	DO I=1,MA
	  WRITE(*,52) I,A(I),SIGA(I)
	  ENDDO

52	FORMAT(I2,4X,F14.6,2X,F14.6)
53	FORMAT(' I         A           SigA')



	CALL MODEL(MA,A,SIGA)



	STOP

	END



